Individual based and mean-field modelling of direct aggregation 
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Abstract. We introduce two models of biological aggregation, based on randomly moving particles with 
individual stochasticity depending on the perceived average population density in their neighbourhood. In 
the first-order model the location of each individual is subject to a density-dependent random walk, while 
in the second-order model the density-dependent random walk acts on the velocity variable, together with 
a density-dependent damping term. The main novelty of our models is that we do not assume any explicit 
aggregative force acting on the individuals; instead, aggregation is obtained exclusively by reducing the 
individual stochasticity in response to higher perceived density. We formally derive the corresponding 
mean-field limits, leading to nonlocal degenerate diffusions. Then, we carry out the mathematical analysis 
of the first-order model, in particular, we prove the existence of weak solutions and show that it allows 
for measure-valued steady states. We also perform linear stability analysis and identify conditions for 
pattern formation. Moreover, we discuss the role of the nonlocality for well-posedness of the first-order 
model. Finally, we present results of numerical simulations for both the first- and second-order model on 
the individual-based and continuum levels of description. 
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1 Introduction 



Animal aggregation is the process of finding a higher density of animals at some place compared to the 
overall mean density. Its formation may be triggered by some environmental heterogeneity that is attrac- 
tive to animals (the aggregate forms around the environmental template), by physical currents that trap 
the organisms through turbulent phenomena, or by social interaction between animals [19] . Aggre- 
gation may serve diverse functions such as reproduction, formation of local microclimates, anti-predator 
behaviour (see for instance [32] for a study of reducing the risk of predation to an individual by aggrega- 
tion in Aphis varians), collective foraging and much more (see, e.g., [14] for a relatively recent survey). 
Aggregation also plays an important role as an evolutionary step towards social organization and collective 
behaviours |15j . These aspects explain the continuing interest in understanding not only the function of 
animal aggregation, but also the underlying mechanisms. The development of quantitative mathematical 
models is an essential part in this quest. While such models first of all help to link individual behavioural 
mechanisms to spatio-temporal group patterns, they also often aim at explaining the observed dynamic 
efficiency of animal groups to adapt to environmental variation, and, moreover, provide a valuable tool 
to study in more detail the system dynamics and their robustness to variations in individual behavioural 
parameters or external parameters, see, e.g., [llj . 

Going back to the pioneering work of Skellam [29], continuum spatiotemporal population dynamics 
are traditionally modelled by reaction-convection-diffusion PDEs and systems thereof. In these models, 
diffusion typically describes the avoidance of crowded areas by the individuals, and as such acts as an 
anti-aggregative force, working against the typically aggregative effect of convection, see the surveys |22j 
and [21J . Our work goes into the opposite direction and shows that biological aggregations can be 
a consequence of solely random, diffusive motion of individuals, who respond to the local population 
density observed in their neighborhood by increasing or decreasing the amplitude of their random motion. 
This kind of behaviour was observed in insects, for instance the pre-social German cockroach Blattella 
germanica. These animals are known to be attracted to dark, warm and humid places |26j . However, the 
works by Jeanson et al. [9j [10] have shown that cockroach larvae also aggregate in the absence of any 
environmental template or heterogeneity. In this case aggregation is the result of social interactions and 
happens as a self-organized process. The individual based model developed and parameterized in [10] 
identified a simple decisive mechanism that can be summarised in the following way: cockroaches do not 
rest for a long time in places with few conspecifics, and once moving they stop preferentially in places 
of high cockroach density. A mathematically better tractable version of this decisive mechanism is the 
model of individuals performing random walks with density dependent coefficients. In the continuum 
limit, this leads to degenerate diffusion models, which, however, depending on their parameters and data, 
might have ill-posed regimes. In particular, Turchin [31] derived a ID model for a population with density 
u{x,t) of the form 



where 4>(u) = (fi/2)u — k$u 2 + (2/co/2u;)u 3 with the positive coefficients \x, ko and oj. Turchin used the 
above equation to describe the aggregative movement of Aphis varians, a herbivore of fireweeds {Epilobium 
angustifolium) . He also discussed the possible ill-posedness of some initial and boundary value problems 
associated with (jl.ip . Depending on the actual profile of </>, Turchin also classified two different types of 
aggregation. Essentially the same model has been derived independently in [T] as a model of cell motility 
with volume filling and cell-to-cell adhesion. The authors showed that the diffusivity can become negative 




(1.1) 
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if the cell adhesion coefficient is sufficiently large and related this to the presence of spatial oscillations 
and development of plateaus in their numerical solutions of the discrete model. Moreover, they used a 
combination of stability analysis of the discrete equations and steady-state analysis of the limiting PDE to 
gain better understanding of the qualitative predictions of the model. Another work studying an equation 
of the type (|1.2p is [23], where existence and uniqueness of solutions of the initial value problem with 
homogeneous Neumann boundary conditions in a bounded domain of R n was shown and some aspects of 
the aggregating behaviour were studied analytically. 

In [25| , following [18] , the authors provide an alternative derivation of (jl.ip for low population density, 
using a biased random walk approach. Then they study the existence of traveling wave solutions for a 
purely negative or zero diffusion equation with a logistic rate of growth g(u), 



du d 
dt dx 



tit \ du 



+ g(u), (1.2) 



discuss the well- and ill-posedness of certain boundary conditions associated with some purely negative 
diffusion equations with logistic-like kinetic part and provide some numerical examples. 

One possibility to overcome the difficulty of the possible ill-posedness of (jl.ip and ([ 1 . 2 j) is to introduce 
a nonlocality, i.e., to substitute the term A<p(u) by A J with 

J(x,t)= / K(x,y)((>(u(y,t))dy 
Jq 

with a nonnegative kernel K(x,y). A particular choice made in [6] is to define K(x,y) as the Green 
function of (I — AA) for a constant A > 0. Equation (|1.2j) becomes then 

^ = A[I-AA]-V(«) + <?(«), 



which is equivalent to 



^ = A U(u) - \g(u) + A^ ) 



This equation was studied in [23] and can be considered as a model of aggregating population with a 
migration rate determined by <fi and total birth and mortality rates characterized by g. In [23] it was 
shown that the aggregating mechanism induced by (j)(u) allows for survival of a species in danger of 
extinction and performed numerical simulations suggesting that, for a particular version of 4>(u), the 
solutions stabilize asymptotically in time to a not necessarily homogeneous stationary solution. Other 
two works going in this direction are [30] , which we discus later (Section I4.4p . and the model of home 
range formation in wolves due to scent marking [2], 

du . / Dnu 

— = A ' 



dt \l+p/a 

dp . , 

— = u(rf + m(p)) - up . 

Here u(x,t) is the population density of wolves, p(x,t) the density of their scent marks and Dq, 7, a and 
\x positive parameters. The increasing function m{p) describes enhanced scent mark rates in the presence 
of existing scent marks. A nonlocal version of the model with m depending on the averaged version of 
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p is also considered. The authors of [2] show that the model produces distinct home ranges; in this case 
the pattern formation results from the positive feedback interaction between the decreased diffusivity of 
wolves in the presence of high scent mark densities and the increased production of new scent marks in 
these locations. 

In our paper we introduce and study two models where formation of aggregates results from random 
fluctuations in the population density and is supported merely by reducing the amplitude of the individual 
random motion in response to higher perceived density. This leads to nonlocal individual-based and PDE 
models, where the nonlocality stems from calculating the perceived density as a weighted average over a 
finite sampling radius. This is usual in modelling biological interactions, see for instance |12 ^ l30 l l20j and 
many more. Our models were inspired by the above mentioned observations of German cockroach [9"1 IIP], 
but does not aim to be a realistic description of their social behaviour. Instead, we consider our work as a 
proof of concept, where the main characteristic of our models is that we do not assume any deterministic 
interaction between the individuals that would actively push them to aggregate. This approach was 
followed for instance in stochastic run-and-tumble models of chemotaxis [27]. Closely related to our work 
is the series of papers [81116} [T7]. However, only the first-order model under specific conditions was studied 
there. The new aspects contributed by our work are the generality of the models (first- and second-order), 
more rigorous derivation of the mean-field limits based on the generalized BBGKY-hiearchy, and rigorous 
well-posedness and stability analysis of the first-order model. 

Our paper is structured as follows: In Section [2] we introduce two stochastic individual based models, 
where every individual is able to sense the average population density in its neighborhood, and respond 
in terms of increased or decreased stochasticity of its movement. In particular, the diffusivity is reduced 
in response to higher perceived density. We present a first-order model, where the location of each 
individual is subject to a density-dependent random walk, and a second-order model, where the density- 
dependent random walk acts on the velocity variable, together with a density-dependent damping term. 
The advantage of the second-order model is that it is possible to introduce a cone of vision, which depends 
on the direction of each individual's movement. In Section [3] we formally derive the corresponding mean- 
field limits, leading to a nonlocal, nonlinear diffusion equation for the first-order model, and a nonlinear 
Fokker-Planck kinetic equation for the second-order model. Moreover, we show that the diffusive limit of 
the kinetic equation leads to the first-order diffusion equation derived before. Then, in Section HI we show 
the existence of weak and measure- valued solutions of the first order mean-field model and perform linear 
stability analysis of the uniform steady states, finding regimes that correspond to the sought-for pattern 
formation (direct aggregation) . Finally, in Section [5] we present the results of numerical simulations of 
our models. We performed Langrangian simulations of the first- and second-order agent-based models 
in a periodic 2D domain, Eulerian simulations of the first-order mean-field model in ID and 2D periodic 
domains, and, finally, of the second-order mean-field model in a spatially ID periodic domain with ID 
velocity. 

2 The stochastic individual based models 

We consider the set of N € N individuals with time-dependent positions X{(t) £ % = 1, . . . ,N, with 
d > 1. Every individual is able to sense the average density of other individuals in its neighborhood, given 



by 




(2.1) 
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where W{x) = w{\x\) with w : M + — > ~R + is a bounded, nonnegative and nonincreasing weight, integrable 
on M. d . A generic example of w is the characteristic function of the interval [0,i?], corresponding to the 
sampling radius R > 0. The average density $j is then simply the fraction of individuals located within 
the distance R from the i-th individual: 

Mt) = jj#{r,\xi-xj\ < R}- 

For the passage to the mean-field limit N — > oo and the analysis of the resulting equations, we will have to 
pose certain smoothnes assumptions on W, which actually exclude the choice of w to be a characteristic 
function. However, from the modelling point of view this is not a concern. 
Let us now introduce our two models: 

• In the first-order model the individual locations are subject to average density-dependent random 
walk, 

dx i = G(# i )dBl i = l,...,N, (2.2) 

where B\ are independent d-dimensional Brownian motions and G : M + — > IR + is a bounded, 
nonnegative and nonincreasing function. For instance, in the numerical simulations of Section [5] we 
use G(s) = exp(— s/3). However, we as well allow for degeneracy, where G(s) = for all s > so > 0. 

• In the second-order model the individuals are described not only by their locations Xi(t) G W d , but 
also by the velocities Vi(t) € M. d , i = 1, . . . , N. The advantage of this description, in contrast to the 
first-order model, is that every individual has a well defined direction of movement. Therefore, the 
weight W in the calculation of the average densities $j can also depend on the relative angle with 
respect to the individual's direction of movement, 

#i(t) = —Y! l W ( x i- x 3> v i)> ( 2 ' 3 ) 

with W(x,v) = w This is important from the modeling point of view, since we can 

define not only the sampling radius R > 0, but also a restricted cone of vision with angle a € (0, 7r]; 
we set then w(s,z) = X[o,fl](s)X[cosa,i](«)- 

The velocity in our model is subject to a density-dependent random walk and a density-dependent 
damping term: 

dxi = Vidt, (2.4) 
d Vl = -H(<&i)vidt + G(#i)dBl. (2.5) 

The function G is as before, while H : IR + — > IR + is a nonnegative and nondecreasing function. 
The damping term —H($i)vi is introduced in order to slow down the agents' motion when they 
approach a crowded place. Obviously, this mechanism is not only necessary to allow for formation 
of aggregates, but also very natural from the modeling point of view. 



5 



3 The mean-field limits 



3.1 The first-order model 

We start with the derivation of the mean- field limit of the first-order model (|2.2p . Unfortunately, the 
standard framework of BBGKY hierarchies cannot be applied, since due to the structure of the model it 
is impossible to obtain a hierarchy where the evolution of the A;-th marginal is expressed in terms of a 
finite number of higher-order marginals. Therefore, we have to apply the recently developed technique of 
introduction of auxiliary variables [3] and their elimination after the mean-field limit passage. 

In particular, under the additional assumption W £ C l (U. d ), we extend the model by introducing 
the average densities given by (|2.ip as a new set of independent variables, governed by the system of 
stochastic differential equations 

Mi = E - x i) ( G ^) dB i ~ G ^j) d#j) • (3.1) 



Let us point out that the random walks Bj, B- in (|3.1[) are correlated with those in (|2.2p . Using the 
Ito formula, we turn to the equivalent formulation of the stochastic system (|2.2|) . (|3.ip in terms of the 
corresponding Fokker-Planck equation 



df 



N 



01 



1 N N A ( 1 

t=i t=i 1 \ j^i 



2 f N 



1 N Pi 



i=l jj^i kj^i 
1 N ff 

+ ^EEEp (V^fe - x k ) ■ VW{ Xl - Xl )G{^ff N ) 

i=l j^i k^i 1 
1 N f& 

i=l jj^i 1 
1 N B 2 

+ ^ E E E QA&w- ^ w{x * ~ Xk) • vw{x i ~ ^wotff") 

i=l jjti kj^i,j 1 3 
1 N r) 2 

■m E E E tm^t - x ^ ■ - ^)G^i?f N ) , 



N2 ^^ t? .Mid#i 

1=1 jjL t k^l,J 



where / = / (t, x\, . . . , xn, $n) is the iV-particle distribution function. Defining the /c-particle 

rN 



marginals by 



fk(t,xi,&i,...,x k ,<& k )= / / / Ar (t,xi,^i,...,XAT,i?Ar)d^ fc+ i...?? A rdx fc+ i... dx N . 

' '(0,oo) JV - fc 



/ / 

jRd(N-k) ^(o 

we obtain the so-called BBGKY-hierarchy for the system (fk)k=v ^ n particular, the equation for f\ 
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flit, x, i?) reads 
dt 2 



1 -& X {G(#?f?) + -^V x • (c(tf ) 2 VW(x - y)fi f (x, 0, y, a) da dy 

+ lw {°^ )2 La L Jo lo VW{X ~ V)VW{X ~ Z)f ^ X ' V ' a ' *' T) dU dT dlJ dZ ) 

+ ^^Uj™G(#) 2 \VW(x-y)\ 2 f»(x,i!>,y,a)dadyY (3.2) 



where we, for the sake of legibility, dropped the indices at x\ and d\. Now, we pass formally to the limit 
N — > oo, assuming that limAr^oo f£* = f k for all k > 1. Moreover, we admit the usual molecular chaos 
assumption about vanishing particle correlations as N — > oo, 



f k (t, xx, ,x k ,-d k ) = JJ/i(t,Xi,i?j) for all k > 2 



i=l 



Then, one obtains from (|3.2p the one-particle equation 

^ = ±A X (G(0) 2 /i) + • {G($?(VW © h)h) + (G(#f\VW © A| 2 /i) , (3.3) 

where the operator © is defined as 

VW®h(t,x)= f I™ VW(x-y)h{t,y,<d)d$dy. 

JR d Jo 

Finally, we reduce (j3.3j) to obtain the standard mean-field description of the system (|2.2p by removing the 
auxiliary variable r &. Indeed, a relatively lengthy, but straight-forward formal calculation shows that 
posesses weak solutions of the form 



/i = g(t,x)S(<& - W* g(t,x)) , with W * g(t, x) = \ W{x-y)g(t,y) 



dy 



and g = g(t, x) satisfies the nonlinear diffusion equation 



^ = ±A X (G(W *g?g) . (3.4) 

3.2 The second-order model 

With the same procedure as before we derive the formal mean-field limit of the model (|2.4p ~ (|2.5p . We 
omit the details here and immediately give the resulting kinetic Fokker-Planck equation for the particle 
distribution function / = f(t,x,v), 



df 
dt 



+ v-V x f = V v - (h{W © f)vf + l -V v (G(W © /) 2 /)) , (3.5) 



where, with a slight abuse of notation, the convolution operator © is defined as 

W®f{t,x,v)= / W(x-y,v)f(t,y,w)dwdy. 
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Let us make the following observation: If the weight W does not depend on v, such that W © / does 
not depend on v as well and W © /(£, x) = W * g(t, x), we can write the (non-closed) system for the mass, 
momentum and energy densities 

g(t,x)=[ f(t,x,v)dv, Qu(t,x)=[ f(t,x,v)vdv, gE(t,x) = - [ f(t, x, v)\v\ 2 dv , (3.6) 

JR d JR d z JR d 



associated with the solution / of (j3.5[) as 



d(gu] 



^ + V x -(gu) = 0, (3.7) 
+ Vx '(/ •^' x '^ lMg> vdv ) = ~ h ( w *8)qu, (3.8) 



dt 

&<kQE) +V x - (±J^f(t,x,v)\v\ 2 vdv\ = -2H(W*g)gE + ^G(W*g) 2 g. (3.9) 



Of 



We observe that only mass is conserved, while momentum and energy are not (neither locally nor globally). 
Indeed, the momentum is dissipated due to the "friction" term, whose strength depends non-locally on g 
due to H(W*g). The energy is, on one hand, dissipated due to the same term, on the other hand is created 
due to the diffusive term in (|3.5p . at the rate ^G(W * g) 2 . In equilibrium, we have H(W * g)gu = 0, 
which means that we either have empty regions (g = 0) or regions with positive density, but zero velocity. 
In these populated regions the equilibrium energy is given by 

F - F ln] - d G{W * Q)2 

3.3 Diffusive limit of the second-order model ( 13. 5 j) 



We show that the first-order model (13.4j) is obtained from (|3.5h in the formal diffusive limit, under the 
assumption that the weight W does not depend on v, which we make for all of this section. Let us recall 
that then W © / does not depend on v as well and W © f(t, x) = W * g(t, x), with g defined by (I3.6j) . 

We start by observing that the equilibria of the collision operator of (|3.5p are given by the local 
Maxwellians 

Therefore, at equilibrium the mean velocity u vanishes if g ^ 0, and the "statistical temperature" T, 
defined by 

Lt=\I M[g](v)\v\ 2 dv, 

z z JR d 

depends non-locally on g and is given by 

Consequently, in crowded regions (where W * g is high) the temperature is low ("freezing of the aggre- 
gates"), while the unhabitated regions are "hot". 
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We introduce the diffusive scaling with the small parameter e > to (|3.5I) . 

e 2 ^ + ev- V x f = V v ■ (h{W * g)vf + l -V v (G(W * g) 2 f)^j , 

and consider the Hilbert expansion in terms of e, / = M[g] + £g, where A4[g] is given by f|3. lOj) . Moreover, 
we perform the Taylor expansion of H{W * g), which we formally write as 

H(W*g) = H + eH 1 + O(e 2 ), 

with 

H = H(W*g M ), H 1 = H'(W*g M )W*g g , 

and similarly for G(W * g) = Go + eG\ + 0{e 2 ). Here, g_M and g g denote the velocity averages of the 
lowest order terms in the expansion, i.e. 



dv. 



Qm= M[g]dv, g g = g 
Then, collecting terms of order e, we obtain 

v ■ v x M = V„ • (H vg + H lV M + ^V,, (G g + G X M) 
which yields, after multiplication by v and integration, 

%V x (qmT\W*qm]) = V«- / Mv®vdv = -H I gvdv, (3.12) 
with T[VK* qm\ given by (|3.1ip . Collecting terms of order e 2 and integrating with respect to v, we obtain 

^tr- + Vx- [ gvdv = , 

Ot JRd 

and using (I3.12|) . we finally obtain the nonlinear diffusion equation 

t V *-(-57^r ?V x (T\W*q M ]qm)) =0. (3.13) 



dt 2 \H(W*g M ) 
Observe that with the choice H = const., f)3. 13j) reduces to (|3.4j) . possibly up to a linear rescaling of time. 



4 Mathematical analysis of the first-order model 

In this section we show the existence of weak solutions of the first-order nonlinear diffusion equation (|3.4p 
and some asymptotic properties of these solutions. For simplicity, we consider the full space setting $7 = 
W 1 , d> 1, but our analysis can be easily adapted to the case of a bounded domain Q with homogeneous 
Neumann boundary conditions, or to the case of periodic boundary conditions, as in the numerical 
examples of Section [5j 
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To simplify the notation, we set F(z) := ^G(z) 2 , thus we rewrite f|3.4j) as 

^ = A(F(W*g)g) , (4.1) 

subject to the initial condition 

e {t = 0) = g . (4.2) 

For the rest of this Section, and without further notice, we make the following reasonable assumptions on 
F and W: 

• F : R + — > R + is a bounded, nonnegative and nonincreasing function. Let us note that we allow for 
degeneracy in (I3.4p . i.e., there might exist an so > such that F(s) = for all s > so- 

• F is continuously differentiable with globally Lipschitz continuous derivative and G = y/2F is 
globally Lipschitz continuous. 

• We W 1,00 (M d ) C\H 



Definition 1 We call g G L°°(0, T; where V{£1) denotes the set of probability measures on $1, a 

weak solution of (|4.ip subject to the initial condition (j4.2j) with 0< a 6 L 2 (Q)r\V(Q), ifV(F(W*g)g) G 
L 2 (0, T; L 2 (f2)) and for every smooth, compactly supported test function ip G C£°([0, T) x fi) u;e /iaue 

I [ g^fdxdt+ [ g ip(t = 0)dx= [ [ V{F(W * g)g) ■ Vipdxdt . (4.3) 
Jo Jn &t Jn Jo Jn 

The proof of existence of weak solutions in the sense of Definition Q] will be performed in two steps. 
First, we consider an approximating, uniformly parabolic equation, and prove the existence of its solutions. 
Then, we remove the approximation in a limiting procedure, to obtain global in time distributional 
solutions. 

4.1 Approximation 

In order to obtain a uniformly positive diffusion coefficient, we use the approximation F £ (z) := F(z)+e for 
e > 0, and analyze the approximating equation S| = A (F £ (W * g)g), which we write in the Fokker-Planck 
form 

^ = V- (gVF £ {W*g) + F £ (W*g)Vg) (4.4) 

subject to the inital condition 

g(t = 0) = f>o • (4.5) 
Theorem 1 For every e > 0, T > and go G L 2 (Q) nV(£l) there exists a nonnegative weak solution 

g £ G L 2 (0, T; H 1 ^)) n ff x (0, T; tf" 1 ^)) n L°°(0, T; 7>(fi)) 
o/ d^.^[ ) in f/ie sense of Definition^ (with F £ in place of F). 
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Proof: The solution is constructed as a fixed point of the map : u i— > g, defined on the convex, bounded 
set 

B R ■= {u G L 2 ((0, T) x fi) ; ||u|| £ o O ( 0i T;£a(n)) < R and G for a - e - * G (°> T )l 
with i? > to be specified later, where £ = 0(-u) is the unique solution of the Fokker-Planck equation 

dg 
di 



V • (gVF £ (W * u) + F £ (W * u)Vg) . 



(4.6) 



Due to the assumed properties of W and F, the diffusion coefficient satisfies F £ (W * u) G L°°((0,T) x $7) 
and is bounded below by e, and the convection coefficient VF £ (W * u) G L°°((0, T) x 0). Thus, (|4.6p is 
a linear, uniformly parabolic equation with bounded coefficients, and the standard theory [5] implies the 
existence and uniqueness of the global, nonnegative weak solution 



g G L°°{0,T;L 2 {n)) D L 2 {0,T; H 1 {n)) n ^(O, T; tf" 1 ^)), 



(4.7) 



such that g(-,t) G V(£l) for almost every t > 0. Therefore, is well defined on Br and, as can be easily 
proven, is continuous there. 

In order to apply the Schauder fixed point theorem (cf. [1]), we need to show that maps Br into a 
compact subset of itself. The key point is the following a-priori estimate, obtained by using g as a test 
function for (|4.6p (note that due to (|4.7p . g is indeed an admissible test function): 



o Jtt 



g^- dx dt + / £>q dx + 
dt 



o Jtt 



F £ \Vg\ 2 dxdt 



o Jtt 



gVF £ ■ Vgdxdt, 



for s G (0, T), where we introduced the shorthand notation F £ := F £ {W * u). Now we use the identity 

gd l dxdt =lljt (1 f dx ) dt =li s)2 dx - \ L & dx 



io Jtt 

and the Young inequality 



o Jtt 



gVF £ - Vgdxdt 



1 

< - 
~ 2 



o Jtt 



F £ \Vg\ 2 dx dt + 2 



I'l 




lo Jtt 





to obtain 



g(x, s) 2 dx + 



1 



o Jtt 



F £ \Vg\ z dxdt< / g^dx 



o Jtt 



VJF F 



g 2 dx dt 



g 2 dx dt . 



Finally, since F £ > e and \F'(y)\ < C for < y < sup a . gC W * u(x) < \\W\\ LO 



V^F £ 
and we conclude 



V^/F £ ~(WWuj 



2 \F'(W*u) 



C 2 

\VW*u\ 2 < — II^H^.oo 



4s 



4F £ (W*u) 

[ g(x,s) 2 dx + ^ [ I \Vg\ 2 dxdt< [ gldx + C £ [ [ g 2 dxdt. 
Jtt 2 Jo J n J n J 

The Gronwall inequality yields an a-priori bound on g in L°°(0,T;L 2 (n), which implies that indeed 
maps Br into itself, for a suitable R > 0. Moreover, reinserting into the right-hand side, we obtain an a- 
priori bound in L 2 (0, T; i^ 1 (SI)) , and an analysis of the right-hand side of (|4.6|) immediately yields a bound 
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ini? 1 (0,T;i7 The proof is concluded by an application of the Aubin-Lions Lemma (cf. [1]), which 

states that the space L 2 (0,T; H 1 (Q)) n H 1 ^, T; .ff" _1 (fi)) is compactly embedded into L 2 ((0,T) x 0), so 
that maps into a compact subset of itself. 



Most a-priori estimates on g £ that have been obtained along the lines of the above proof are not 
uniform in e and can thus not be used to pass to the limit e — > 0. However, we can tackle this problem 
with a more careful analysis, where the key idea is to look for estimates on ^/F^g 6 , rather than on g £ . 

Lemma 1 There exists a constant C independent of e > small enough, such that the solutions g £ of 
\4-4\)> constructed in Theorem^ satisfy 



l0 £ |ll,°°(O,T;-p(fi)) - 1j \\Q £ \\l° c {0,T;L 2 (Q:)) — @ J 



dg £ 



Of 



<c. 



L 2 {0,T;H-1{£1)) 



and 



\W * Q e \\L^{o,T;W^°°(n)) - C ' 



w * 



dg £ 
~dt 



< C . 



L°°(0,T;L 2 (fl)) 



Proof: For the sake of better legibility, we will drop the superscript at g £ in the proof. Since g is 
constructed as a nonnegative weak solution of the Fokker-Planck equation (|4.4[) , the first estimate follows 
immediately due to the mass conservation 



lle(V)ILi(n) = / Q(t,x)dx= / g (x)dx = l 
Jn Jn 



Proceeding in a similar way as in the proof of Theorem [H we conclude for s £ (0, T) 

\ [ g(s,x) 2 dx+ [ [ F £ \Vg\ 2 dxdt + [ [ gVF £ • Vgdx dt = [ g 2 dx , 
2 Jn Jo Jn Jo Jn Jn 

where we again use the shorthand notation F £ := F e (W * g). Using the identity 



V y/F e g 



e^VFe 



gVF £ -Vg + F £ \Vg\ 2 

we obtain 

\ [ g(s,x) 2 dx+ [ [ v(^f¥ £ g 
1 Jn Jo Jn v 

Now, since \G'(y)\ < C for < y < sup^gQ W * u(x) < \\W\\ LOO , we have 



dxdt =}- ( Qq dx + l 

2 Jn Jo Jn 



VWF £ 



dx dt . 



V^/F £ (W*g) = \G'(W*g)\\VW*g\ <C\\W\\ w i l00 . 



Hence, 



l - [ g(s,x) 2 dx+ ( [ V(^/F £ g) 2 dxdt<l [ g 2 dx + C 2 \\W\\ 2 wl , x [ [ g 2 dx 
* Jn Jo Jn v ' ^ Jn Jo Jn 



dt. 
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and the Gronwall inequality yields a uniform estimate for g in L°°(0, T; L 2 (Q)), which subsequently implies 
a uniform estimate for \J F £ (W * g)g in L 2 (0, T; H 1 ^)) and, subsequently, also for ^| in L 2 (0, T; 



Finally, we derive the estimates for W * g. Since W G VF 1,00 (lR d ), we immediately obtain 

II w * e\\ L °°(o,T-w^°°(p.)) ^ II^IIm/i.-^^^sup ^ dx = ||Ty|| w i j00(Rd) . 

Moreover, we have 

w * Tt = w * A(i?£ ^ = AW * (Fe£,) ' 

and with W G i? 2 (M d ), G L°°(0, T; P(O)) and the uniform boundedness of F £ in L°°((0,T) x Q), we 
obtain a uniform bound for W * (|f ) in L°°(0, T; L 2 (0)). 

■ 

4.2 Global existence 

From the above approximation and a-priori estimates we can easily pass to the limit e — > and conclude 
the existence of a weak solution for (|4.ip : 

Theorem 2 Let go G L 2 (f2) nT'(fi) and T > 0. T/ien i/iere exists a solution 

g G L°°(0, T; L 2 (0)) n L°°(0, T; n i/^O, T; iT 1 ^)) 

°f (|4.ip ~ (|4.2p in i/ie sense of the weak formulation (|4.3j) . suc/i i/iai in addition 



^F(W*g)geL 2 {0,T;H 1 (Q)). 

Proof: We use the uniform estimates of Lemma Q] and the Banach-Alaoglu theorem [3] to extract a 
subsequence g £n such that 

g £n ->>* g weakly-* in L°°(0, T; L 2 (0)) , 

^ ^ |f weakly in L 2 (0, T; F" 1 ^)) , 

y/F^W * QF")^" n weakly in L 2 (0, T; ff 1 ^)) , 

for some u G L 2 (0, T; i7 1 (r2)). Moroever, we use a variant of the Aubin-Lions Lemma [28] to conclude 
the compact embedding of L°°(0,T; W 1,0 °(fi)) n M /1,oo (0, T; L 2 (0)) into L°°(0, T; (7(0)). Thus, we can 
extract a further subsequence, again denoted by g £n , such that as e n — > 0, 

W * g £n -> v strongly in L°°(0, T; (7(0)) . 

Since further 

W * £ £n ^* W * g weakly-* in L°° (0, T; L 2 (fi) ) , 
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we conclude v = W * g by the uniqueness of the limit. Due to the continuity properties of F and F' , we 
also have 

F £n {W * g £n ) -> F{W * q) strongly in L°°(0,T; C{Q)) , 



V F e n (W*Q^) -> v 7 ^^*^) strongly in L°°(0,T; C(fi)) , 

F e ' n (W * e e ") -»• * strongly in L°°(0, T; . 



Consequently, 



y/F £n (W*f")g Sn ->> v / F(W*g)g weakly-* in L°°(0, T; L 2 (0)) 



and, again, by the uniqueness of the limit we identify u = y/F(W * g)g. Finally, we use the identity 



V (F En (W * g £ ") g e ") = F' £n (W * g £ ") g £ " V (W * g £n ) + V ' F £n (W*g^)V (V F £n (W * ^) g £n 
and the above limits to conclude 

V (F £n (W * g £n ) g £n )^V (F {W * g) g) weakly in L 2 (0, T; L 2 (tt)) . 



Hence, we can pass to the limit e n — > in the weak formulation (j4.3j) and conclude that g is a weak 
solution of (IQl-dOl). 



Finally, we want to reduce the regularity of the initial value towards solely probability measures. This 
is relevant since we shall see below that indeed Dirac 5 distributions can be stationary solutions. For this 
sake we define the very weak (distributional) notion of the solution: 

Definition 2 We call g G L°°(0,T;V(CI)) a very weak (distributional) solution of (|4.ip subject to the 
initial condition (I4.2|) with go G V{0), if for every smooth, compactly supported test function (p G 
C c °°([0,T) x fi) we have 

[ [ *^fgdxdt + [ <p(t = 0)g dx = - [ [ {F{W * g)Aip)gdx dt , (4.8) 
Jo Jn Jq J Jq 

where we denote by gdx the integration with respect to the probability measure g(t,-). 
Theorem 3 Let go G V(£l) and T > 0. Then there exists a very weak solution 

geL°°(0,T;V(n)) 

of P~T1) - P~2T) in the sense of (|Q]> . 

Proof: Let us consider a sequence (£»o) n gN — L 2 {Q) n V(Q) such that Qq — > go in V(£l). Moreover, let 
us denote by g n the corresponding weak solutions of (|4.ip ~ (|4.2p . which are trivially also the very weak 
solutions in the sense of (|4.8p . Conservation of mass provides immediately a uniform bound on g n in 
L°°(Q,T;V(Q)), such that for a subsequence, again denoted by g n , we have 

g n g weakly-* in L°°(0, T; V(£l)) . 
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To show that g is a solution to (|4.8p . we need to prove that F(W * g n ) converges to F(W * g) strongly 
in L 1 (0, T; C(fi)). First, let us note that, due to the assumptions on W, W * g n is uniformly bounded 
in L°°(0, T; W 1,oc (fl)), and the same holds for F(W * g n ). Consequently, the sequence F(W * g n )g n is 
uniformly bounded in L°°(0,T;V(Q)). From (jHJ) it immediately follows that ^ = A(F(W * is 
uniformly bounded in L°°(0, T; (C% (fi))*), where (C^(fi))* denotes the dual space to C^(O), the space 
of twice continuously differentiable functions with compact support in 0. Consequently, the sequence 
jf t F{W * g n ) = F'(W * g n )W * ^ is uniformly bounded in L°°(0,T; (C c 2 (0))*), and the generalization 
of the Aubin-Lions lemma by Simon [28] implies then the strong convergence of F(W * g n ) to F(W * g) 
in L 1 (0,T;C(fi)). 



4.3 Stationary Solutions 

The numerical simulations provided in Section [5] suggest that for G > (for instance, G(s) = e~ s ) and 
bounded domain with periodic boundary conditions, the stationary solutions of (|3.4p consist of one or 
more localized, but not compactly supported aggregates (clumps), see the bottom right panel of Fig. [5] 
and Fig. [7J However, we were not able to characterize these stationary aggregates analytically. Instead, 
we provide a few examples of other types of stationary solutions, posed either in the full space setting 
12 = R d or on a torus Q = T d with periodically extended W. These examples are rather trivial, however, 
still provide an interesting insight into the relatively rich structure of the solutions of (|3,4p . 

• The most trivial type of stationary solution is the constant state g = c > 0. Clearly, this has finite 
mass only if $7 = T d . 

• If there exists an sq > such that G(s) = for all s > sq, then any profile g such that 
Q ^ s (/o W(x) dx) almost everywhere on is a stationary solution to the distributional formu- 
lation (]4.8p . Indeed, such a solution satisfies W * g > Sq, so that GiW * g) = 0. However, again, 
this solution has infinite mass if O = M. d . 

• If there exists an so > such that G(s) = for all s > sq, and, moreover, W is continuous on Q, 
we construct the atomic measure 

N 

Q( x ) = y^ y (HS(x - Xi) 
1=1 

for some Af £ N, £ and q > such that CiW(0) > so for all i = 1, . . . , N. Then g is a 
distributional stationary solution to (14. 8p . Indeed, for any i = 1, . . . , N we have 

N 

W * g( Xi ) = CjW(xi - Xj) > Ci W(0) > s . 
i=i 

By the continuity of W and, consequently, W * g, we have W * g > so on some neighbourhood 
of Xi. Therefore, G{W * g) = on some open set containing U i= i x i an di finally, G(W * g)g = 
everywhere on 12. 
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4.4 Linear stability analysis 

In this section we perform a linear stability analysis of constant density states of the nonlinear diffusion 
equation (|3.4p . As discussed previously, we have constant steady state solutions either in the full-space 
setting = M. d (however, with infinite mass), or on a bounded domain $7 C K with periodic boundary 
conditions. Without loss of generality, we assume the normalization J^W(x) dx = 1. Let us make the 
perturbation ansatz g = go + eg, where go > is a constant state such that G(qq) > and G'(go) < 0. 
Then W * g = go + eW * g and assuming sufficient smoothness of G, we have 

G(W * q) 2 = G(q ) 2 + 2eG(Q )G'(g )W * g + 0(e 2 ) . 
Inserting the ansatz into (13.4)) and collecting terms of order 0(e), we arrive at 

dg_ _ G|o) (G( ^ o)A ^ + 2G>(q )q A(W *q))=0. 
Performing the Fourier transform, we obtain 

dQ_ + (o {go) + 2G'(g )g w) g = 0, (4.9) 

where we denoted g = g(t,£) the Fourier transform of g. Consequently, with the assumption G(go) > 
and G'(go) < 0, those wavenumbers £ of g are stable for which 

2G'{g )go 

Since W E L 1 (M d ), we have W S C°(]R d ), and, therefore, all wavenumbers of g larger than a certain 
threshold will be stable. On the other hand, wavenumbers violating (|4.10p will lead to pattern formation, 
as we show in our numerical examples, Section [5j Moreover, a quick inspection of (|4.9p leads to the 
expectation that the stable high wavenumbers will be smoothed-out on a faster time scale, compared to 
the slower emergence of patterns due to the unstable lower wavenumbers. Finally, considering the scalings 
W r (x) := r~ d W\(x/r) of a fixed kernel W\ G L 1 (]R' i ), we have W r {£) = W\(r£) and, therefore, we expect 
that the wavelength of the patterns (size of aggregates) will scale with r. This can be as well clearly 
observed in our numerical examples. 

It is interesting to consider the formal extremal case W = 5q, i.e., W * g = g. Then W = 1 and we 
conclude that the constant state £o is stable if and only if 

(G^o) 2 ^)' = G(g ) 2 + 2G(eo)G' / (£o)£o > • ( 4 - u ) 

In fact, violation of (14.111) with W = 5q leads to an ill-posed problem, since then (|3.4j) looks like a backward 
heat equation, which can be seen if we expand the derivatives and write it in the Fokker-Planck form as 

|f = |V- (G(g)[2G'(g)g + G(g)]Vg) . (4.12) 

The equation is parabolic (and thus well-posed) only if the diffusivity 2G(g)G'(g)g + G 2 (g) is strictly 
positive, which is our stability condition (14. lip . Therefore, only imposing an initial condition go uniformly 
satisfying (|4,lip leads to a well-posed diffusion equation for all t > 0, since (|4.1ip is preserved due to 
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the maximum principle. On the other hand, if G > 0, the nonlocality W G L 1 (R rf ) always stabilizes the 
equation in the sense of (I4.10p . Indeed, writing it in the Fokker-Planck form 

— = —V • iylndividualbasedandmean — fieldmodellingofdirectaggregation [2G(W * g)G'(W * g)\7W * g\ g + G 2 (\ 

we see that the second-order term appears with the positive diffusivity G 2 (W*g). The first-order term de- 
scribes then the transport of g along the generalized gradient VW*g and is responsible for the aggregative 
effect. 

Clearly, the ill-posedness of (|4.12p can be avoided by merely introducing a nonlocality in the first-order 
transport term, while the diffusivity may stay local (and possibly degenerate). Such a model was derived 
and studied in [30], which with our notation is written as 

^ = V-(-eVW*e + e 2 V e ) . 

This equation was constructed as a model for biological aggregations in which individuals experience 
long-range social attraction and short range dispersal. Let us note that here, in contrast to our model, 
the diffusivity is an increasing function of the density g. In [30] it was shown that it produces strongly 
nonlinear states with compact support and steep edges that correspond to localized biological aggregations, 
or clumps. Similarly as can be observed in our numerical simulations in Section these clumps are 
approached through a dynamic coarsening process. 

Another insight into the stabilizing effect of the nonlocality is provided by the introduction of a formal 
expansion of W. Taking W as the standard mollifier and W e (s) := e~ d W(s/e), such that W £ — > 5q as 
e — > 0, we expand 

W £ *g = e~ d [ w(^l) e (y)dy 

JWL d V 6 / 

= / W(-z)g(x + ez) dz 
Jm. d 

= f W(-z)( g(x) + eVg(x)z + -e 2 z t V 2 g(x)z) dz + 0(e 3 ). 

jR d V 2 / 

Now, due to the symmetry W(—z) = W(z) and the normalization f Rd W(z) dz = 1, we have 

W £ *g = g+ £ -^/3Ag + 0(e i ), 
where the constant j3 > is such that f Rd W(z)ziZj dz = (3Sij. Inserting this into (|3.4|) . we obtain 

dtQ = ^A^G^g+ £ ^Ag + 0(e 4 ^ 

= (G(g) 2 g + e 2 f3G(g)G'(g)gAg) + 0(e 4 ) . 

Up to the terms of fourth order in e, this is a Cahn-Hilliard-type equation which is well-posed for every 
e > since the term e 2 f3G(g)G' (g)g is strictly negative if G(g) > and G'(g) < 0. However, if e = 0, i.e. 
W = 5q, this regularizing effect is lost. 
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Figure 1: The first-order individual based model with N = 400 agents and sampling radius R = 0.025, 
subject to a random initial condition. 



5 Numerical examples 

In this section we present numerical examples for the first and second-order individual based models (12.2p , 
(|2.4p - (|2.5p in 2D, the first-order mean-field limit (|3.4j) in ID and 2D and the second-order mean-field 
limit (|3.5p in ID with ID velocity. 



5.1 The first-order individual based model ( 12.21) 

We consider a system consisting of N = 400 individuals in a 2D domain £1 = (0, 1) x (0, 1) with periodic 
boundary conditions. The initial positions are generated randomly and independently for each individual 
from the uniform distribution on we took the same initial condition for all the three experiments below. 
We choose W(x) = w(\x\) with w the characteristic function of the interval [0, R], corresponding to the 
sampling radius R, and for R we take the values 0.025 (Fig. [[]), 0.05 (Fig. [2]) and 0.1 (Fig. [3]). For G 
we make the choice G(s) = exp(— s/3). The system of stochastic differential equations ()2.2[) is integrated 
in time using the Euler-Mayurama scheme with time-step length t = 10 -3 . We used the linear stability 
analysis in Section [5] to make the "right" choice of G, such that we could observe pattern formation. 
Indeed, if G decreases too quickly, the system will "freeze" immediately, before any aggregates could be 
formed; on the other hand, if G does not decrease fast enough, the system is "overheated" and does not 
allow aggregates to persist. 

In Figs. Q] [2] and [3] we observe that with a smaller sampling radius R, a larger number of small 
aggregates is created on a faster time scale. The choice R = 0.1 (Figj3]) leads to creation of one sigle 
aggregate. This aggregate is approximately ring-shaped, with higher density of particles around the 
circumference and lower density in the middle. This can be explained by the fact that the aggregate 
grows by "capturing" particles from its neighborhood, and once a particle is captured, its mobility is 
greatly reduced, so that it only slowly makes its way towards the center of the aggregate. Let us also 
mention that the right-most panels in Figs. [THS] present quasi-steady states, where the aggregates are 
in a dynamic equilibrium with a very few free-running particles. However, on a very long time scale, 
the smaller aggregates typically disintegrate and their particles are caught by the larger ones. Such a 
coarsening behavior is typical to diffusive aggregation systems, as for instance the Cahn-Hilliard equation, 
or the nonlocal continuum model of }30] , 
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t=0.1 t=0.25 t=0.7 



Figure 2: The first-order individual based model with N = 400 agents and sampling radius R = 0.05, 
subject to a random initial condition. 




Figure 3: The first-order individual based model with N = 400 agents and sampling radius R = 0.1, 
subject to a random initial condition. 



5.2 The second-order individual based model (12.4ft — ( 12751 ) 

Generally, the behavior of the second-order individual based model is very similar to this of the first-order 
one, at least if we do not impose a restricted cone of vision (i.e., W in (|2.3p does not depend on v). 
Indeed, with a suitable choice of parameters, one again observes the formation of quasi-stable aggregates, 
whose number and size depend on the sampling radius. The most striking difference with respect to the 
first-order model is that the movement of the agents is smoother (their velocities are continuous) and the 
shape of the aggregates is more complex (we observed emergence of ellipsoidal aggregates, instead of the 
almost-circular ones in the first-order model). 

The situation becomes slightly more interesting if we consider a restricted cone of vision - we choose 
180°, such that W(x, v) = w (\x\, with w(s, z) = X[o,R]( s )X[o,i]i z ) ■ The sampling radius is chosen as 

R = 0.05. Again, we simulate with N = 400 individuals in a 2D domain O = (0, 1) x (0, 1) with periodic 
boundary conditions, and G(s) = exp(— s/3). For H we set the constant H = 2. The system of stochastic 
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Figure 4: The second-order individual based model with N = 400 agents, sampling radius R = 0.05 and 
180° cone of vision, subject to a random initial condition. 



differential equations (|2.4f ) — (j2.5f) is integrated in time using the Euler-Mayurama scheme with time-step 
length t = 10~ 3 . The initial positions of the agents are generated randomly in the same way as for the 
first-order model, while their initial velocities are generated independently and randomly from the 2D 
normalized Gaussian distribution. Three snapshots of the evolution of the system are shown in Fig. [U 
the velocities being marked by the linear elements for each individual. The right-most panel in Fig. [3] 
shows the quasi-steady state with two aggregates formed, in dynamic equilibrium with the "free-running" 
individuals. It is obvious that, compared to the previous simulations, the aggregates are less densely 
packed and their shapes are less circular. Also, the portion of the "free-running" particles is much higher, 
which clearly is an effect of the restricted cone of vision (the captured particles can leave the aggregate 
more easily). 



5.3 The first-order mean-field model (13.41) in ID 

We simulated the first-order mean-field model (I3.4p in the ID periodic domain (0, 1), using semi-implicit 
finite difference discretization for the space variable and first-order forward Euler method for the time 
variable. The space grid consisted of 200 equidistant points, the time step was 10 -4 . As before, we chose 
W(x) = w(\x\) with w the characteristic function of the interval [0,0.1], and G(s) = exp(— s/3). We 
imposed a random initial condition for g, generated such that for every grid point a random number from 
the uniform distribution in [0, 1] has been drawn. Snapshots of the evolution are shown in Fig. [5j We 
observe that quite a strong smoothing effect takes place on the fast time-scale (first row in Fig. [5|), while 
aggregation takes place on a time-scale approximately one order of magnitude slower (second row); this 
is explained with the stability analysis in Section I4.4i First, two aggregates of different sizes are created, 
however, both of them are unstable - the smaller one is smoothed out, while the larger grows further, 
until the steady state is reached (lower right panel). Observe also the characteristic "fork"-like shape of 
the profile in the lower mid panel (t = 2.65), which is due to the mass arriving from the neighbourhood 
with higher diffusivity than the diffusivity in the middle of the profile; compare also with the ring-shaped 
aggregate in the right panel of Fig. [2j However, this fork-like structure is eventually also smoothed out, 
to finally obtain the steady profile. Let us note that the steady aggregate, although well localized, is not 
compactly supported, i.e., the profile has positive density everywhere on [0, 1]. 



20 




t=0.8500 t=2.6500 t=11.0500 



Figure 5: The first-order mean-field model in a periodic ID setting with random initial condition (upper 
left panel). Steady state was reached at t ~ 11 (lower right panel). 

5.4 The first-order mean-field model (13.4ft in 2D 

We simulated the first-order mean-field model (|3.4p in the 2D periodic domain (0,1) 2 , using the same 
type of discretization as in the ID case. The space grid consisted of 100 x 100 equidistant points, the time 
step was 10 -4 . We chose the sampling radius R = 0.07, i.e., W(x) = X[o,o.7] (MX an d G(s) = exp(— s/3). 
Snapshots of the evolution are shown in Fig. [6l Starting again from a random initial condition, we 
observed a rapid formation of approximately ring-shaped pre-aggregates, which eventually turn into an 
almost regular pattern of well localized (but not compactly supported) clumps. However, the smaller 
aggregates may be unstable and diffusively disintegrate, and their mass is absorbed by their neighbors, 
as shown on the lower right panel of Fig. [6l Finally, in Fig. [7] we present examples of patterns produced 
with the sampling radii R = 0.6 (left panel) and R = 0.11 (right panel). 



5.5 The second-order mean-field model (13.5ft 

We conclude this section with simulations of the second-order mean-field model (|3.5p in ID. The spatial 
domain fi = [a, b] is discretized using an equidistant mesh x% = a + iAx, the velocity domain V = 
[v -m m , Bmrn] at grid points Vj = v m i n + j Av . The time step At satisfies the CFL condition for f max , i.e., 
At = , The numerical scheme is based on a splitting method: Given an initial datum f(x,v,0) = 

fo(x,v), we split the system at every time tk = kAt into the following steps: 
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0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 



t=0.900 t=0.950 

Figure 6: The first-order mean-field model in a periodic 2D setting with random initial condition (upper 
left panel). After the initial rapid smoothing of the high-frequency components, several ring-shaped 
structures are created (upper right panel), which eventually turn into an almost regular pattern of well 
localized aggregates (lower left panel). However, the smaller aggregates may be unstable and diffusively 
disintegrate, and their mass is absorbed by their neighbors, (lower right panel). 
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Figure 7: An example of patterns produced by the first-order mean-field model in a periodic 2D setting 
with the sampling radii R = 0.6 (left panel) and R = 0.11 (right panel). 



1. Solve transport equation in x 



df 

at 



+ Vi ■ V x f = 0, 



(5.1) 



subject to periodic boundary conditions on $7, for every V{ on the time interval t G [£&,£& + |A£], 
using an upwind scheme with the superbee flux limiter. 



(5.2) 



2. Starting with the solution of the transport equation (|5.ip . solve 

^ = V« ■ U(W ® f)vf + \v v {G(W ® fff) 



with no flux boundary conditions on V, on the time interval [ifc,£fc+i], using a semi-implicit time 
discretization. 

3. Finally solve (|5.ip using the solution of (|5.2p for another half time step |At. 

We set the computational domain to fl = [0, 1], the velocity domain to V = [—1,1] and the mesh sizes 
to Ax = 10 -2 and Av = 2 x 10 -2 . We choose similar conditions as in the individual based model, i.e. 
a limited cone of vision W(x,v) = w(\x\, with w(s,z) = X[0.R) ( s )X[o,i] ( z )- The sampling radius is 

set to R = 0.07, G(s) = exp(— 2s) and H = 2. The initial datum corresponds to a small perturbation 
of a uniform distribution with mass one. Snapshots of the evolution of the particle distribution density 
f(t, x, v) and the mass density p(t, x) at different times are depicted in Fig. We observe a fast smoothing 
of f(t,x,v) in time and a subsequent formation of a stable aggregate. 

If we decrease the sampling radius to R = 0.03, two separate aggregates form, see Fig. [9] Again, we 
observe that with a smaller sampling radius the aggregation happens on a faster time scale. 
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(a) Particle distribution function f(t,x,v) at t = 4 
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(b) Mass density p(t,x) at t = 4 
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(c) Particle distribution function f(t,x,v) at t = 8 (d) Mass density p(t,x) at t = 8 
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(e) Particle distribution function f(t,x,v) at t = 20 (f) Mass density p(t,x) at t = 20 
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Figure 8: Second order model with limited vision and R = 0.07; the red line in the left column indicates 
the initial datum. 
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(a) Particle distribution function f(t,x,v) at t = 4 
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(c) Particle distribution function f(t,x,v) at t = 8 
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(e) Particle distribution function f(t,x,v) at t = 20 
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Figure 9: Second order model with limited vision and 
the initial datum. 
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(b) Mass density p(t,x) at t = 4 
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(d) Mass density p(t,x) at t = 8 
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(f) Mass density p(t,x) at t = 20 
R = 0.03; the red line in the left column indicates 
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